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Abstract. The problem of successfully simulating ionic fluids at low temperature 
and low density states is well known in the simulation literature: using conventional 
methods, the system is not able to equilibrate rapidly due to the presence of strongly 
associated cation-anion pairs. In this manuscript we present a numerical method for 
speeding up computer simulations of the restricted primitive model (RPM) at low 
temperatures (around the critical temperature) and at very low densities (down to 
10 -10 er~ 3 , where a is the ion diameter). Experimentally, this regime corresponds 
to typical concentrations of electrolytes in nonaqueous solvents. As far as we are 
aware, this is the first time that the RPM has been equilibrated at such extremely 
low concentrations. More generally, this method could be used to equilibrate other 
systems that form aggregates at low concentrations. 



PACS numbers: 61.20.Ja, 61.20. Qg 



1. Introduction 



Computer simulations have yielded invaluable insights on the properties of ionic fluids. 
The nature of fluid-fluid ('vapour-liquid') phase separation and the universality class of 
the associated critical point have attracted particular attention. In these studies, the 
restricted primitive model (RPM) has played a central role [1-17]. The RPM is a simple 
representation of molten salts and ionic solutions. It consists of an equimolar binary 
mixture of positively and negatively charged hard spheres with charges ±q and equal 
diameters a, immersed in a continuum with dielectric constant e. In terms of the model 
parameters, the reduced temperature is defined as T* = kBTDa/q 2 , where ks is the 
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Boltzmann's constant, T is the absolute temperature, D = Airee^, and eo is the vacuum 
dielectric permittivity; the reduced density is defined as p* = pa 3 , where p = N/V is the 
total number of ions per unit volume. The most recent high-precision Monte Carlo (MC) 
simulations locate the critical point at a critical temperature T* ~ 0.05 and a critical 
density p* ~ 0.08 [18,19]; the critical point has been confirmed as belonging to the 
three-dimensional Ising universality class [19]. Interestingly, the critical temperature is 
close to room-temperature conditions for sub-nanometre monovalent ions in oily solvents 
with e ~ 5-10. However, the ion concentrations in these nonaqueous electrolyte solutions 
are often in the nM regime (p* ~ 10~ 10 ) [20-22] which motivates the parameter choice 
of the present study. 

It has been clear for a long time that physical clustering of the ions has an 
important effect in the vapour region [4, 23], as strongly suggested by the Bjerrum theory 
[24]. Analysing the features of this system in the low density-low temperature regime, 
Valleau [25] and Gillan [26] showed that the ionic fluid tends to form dimers, triplets, 
and higher order clusters, and that clustering has a crucial effect on the equilibrium 
properties of the RPM. Weis and Caillol [27] and Bresme et al. [28] characterised the 
cluster structures quantitatively at temperatures around T c and at densities around 
p c /5 and p c /3, respectively. Later on, Camp and Patey identified different regimes of 
ion association well below p c [29]: at low temperature the system apparently consists 
of only clusters; at intermediate temperature the system is predominantly associated, 
but with some free ions; and at high temperature the majority of ions are free. The 
RPM with screened Coulombic interactions can serve as a model for charged colloidal 
systems [30-32] . Caballero and coworkers studied such a model with an inverse screening 
length of k = Qa~ l , mimicking the effect of added electrolytes present in the medium: 
the critical point was located at T* ~ 0.17 and p* ~ 0.22, and the familiar clustering 
phenomenon (of the colloids) in the dilute phase was observed. 

The main obstacle to simulating ionic fluids successfully in the low-temperature 
regime, where coexistence occurs, is the strong association of ions at distances close to 
contact and the resulting extremely slow equilibration. Graham and Valleau [8] pointed 
out that, when studying the low-temperature regime, conventional MC or molecular 
dynamics methods are not sufficient to equilibrate the system. Therefore, they first 
used a type of umbrella sampling named "temperature scaling Monte Carlo" at several 
densities [33]. Next, Valleau proposed "density scaling Monte Carlo", a novel algorithm 
based on umbrella sampling over broad ranges of densities [34], and applied it to the 
RPM near the critical point. Orkoulas and Panagiotopoulos computed the vapour-liquid 
phase diagram [11], and to accelerate convergence, proposed ion-pair and cluster moves 
capable of grouping and moving clustered ions. The primary motivation in this work 
was the computation of the vapour-liquid phase diagram, and hence reduced densities 
of no less than 10~ 4 were considered. 

In recent work, Allahyarov and co-workers [35] have shown that, around the 
critical temperature, oppositely charged micro- ions tend to form 'Bjerrum pairs', in 
which oppositely charged particles are closer than the Bjerrum length \ B = q 2 /ksTD. 
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The authors of reference [35] found that at the lowest salt concentrations (around 
10 -9 mol L _1 , corresponding to a reduced density of around 10~ 10 for ions with a — 5 A), 
almost 90% of the ions resided in pairs. In a later publication [36], the same authors 
used a different definition of a cluster (with a cut-off of A = 3er), and narrowed their 
results down to a smaller concentration range: the new results became valid only for salt 
concentrations between 10 -4 mol L -1 and 10~ 2 mol L" 1 (or reduced densities between 
10 -5 and 10~ 3 ). The main problem found by the authors was equilibrating the system 
at extremely low densities by means of standard simulation techniques. 

The aim of our work is to present a novel MC technique that achieves rapid 
equilibration of the RPM at low temperatures (around T c ) and very low reduced densities 
(from 10~ 3 down to 10~ 10 ). This method might also be applied to the equilibration of 
other systems that form aggregates at low concentrations. 



2. Simulations 



The interaction potential for the RPM is 

U[n ' ] = { m,IDr r!, > a, « 
where qi = ±q. The system is comprised of N/2 cations and N/2 anions in a cubic box 
of length L, with periodic boundary conditions (PBCs) applied. We use MC simulations 
of iV = 256 ions in the NVT ensemble. The choice of such a relatively small number 
of particles is justified by the fact that, as the densities under study are so low, the 
simulation box is always large enough to exclude any significant finite-size effects due to 
the PBCs; with N = 256 ions and p* oc 10 -10 , the simulation box length is around 10 4 <r. 
We checked for finite-size effects by running simulations at p* = 10~ 4 with either 256 
or 1000 particles, and making sure that the computed energy per particle was the same 
within statistical uncertainties. Moreover, the box lengths are large compared to the 



range of Debye-like screening, equal to the Debye length \d/o~ = yT* /Anp*. Table 1 
shows the values of L and A^ at the densities and temperatures considered in our work; 
the density range is 10~ 10 < p* < 10~ 3 and the temperature range is 0.04 < T* < 0.07. 
In this paper we will concentrate on simulations at the lowest density and temperature; 
the full range of state points will be considered in a forthcoming publication. 

The long-range interactions were handled using the Ewald sum with tin-foil 
boundary conditions [37-40] . For each density we carefully tuned the Ewald parameters 
a, r c and k max , being the width of the Gaussian distribution characterising the screening 
term in real space, the real-space cut-off, and the reciprocal-space cut-off, respectively. 
a was chosen using the empirical rule aL = 5.6 [41], r c was set to L/2, and /c max such 
that the relative error in the reciprocal-space sum was of the order of 10~ 5 [38]. The 
values of a are indicated in table 1; /c max was always set to 10 x (2-k/ L). To test our 
code, we computed the energy per ion pair u v in a liquid at T* = 0.042 and p* = 0.17; 
we obtained (/3u p ) = —1.26 ± 0.01, which is in perfect agreement with that computed 
for the same state point by Romero-Enrique et al. [42]. We have also computed the 
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Table 1. Box edge L, Debye screening length Ad, and Ewald real-space screening 
parameter a for all of the simulated densities and temperatures ranging from T* = 0.04 
to T* = 0.07, and with N — 256. For each density, the smallest value of Ad corresponds 
to the lowest temperature, and the largest Ad to the highest temperature. 
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energy per particles at rho=0.175 and T=0.05, and compared with the results in Table 
VII of Ref[27]: our results is (U /NksT) = —12.38 ± 0.01 in perfect agreeement with 
their results of U/NksT = —12.38. Moreover, we computed the energy per particle at 
lower densities where the system is in a vapour phase, at p* = 0.002 and T* = 0.05[27], 
and found (U /NksT) = —10.20 ± 0.05, in good agreement with that computed for the 
same state point by Caillol and Weis[27] (U/Nk B T = -10.15). 

It is well known that the RPM forms clusters in the subcritical vapour phase. In 
order to identify the clusters, we use Gillan's definition, according to which two particles 
belong to the same cluster if they are separated by a distance shorter than a given cut- 
off A [26] . In this way, we detect the total number of isolated ions, the total number of 
associated ions, and the total number of clusters of a given size. In what follows, and 
unless stated otherwise, we will study the cluster formation when A = 2a. 

3. Techniques to equilibrate the RPM at low temperatures and densities 

In order to study the RPM at low temperatures and densities, ad-hoc simulation methods 
have been employed to overcome the problem of slow convergence towards equilibrium. 
The main obstacle to simulate the RPM in the low-temperature region is the strong 
binding effect of oppositely charged ions at short distances, as the thermal energy 
available to drive two oppositely charged particles away from each other is much less than 
the attractive Coulomb energy, i.e., T* C 1. Therefore, in order to reach equilibrium, 
the system would have to be simulated for a prohibitively long time. In our simulations, 
this equilibration problem is going to be even more pronounced, since we aim to study 
very dilute systems where isolated ions are so far apart from each other that they spend 
most of the time freely diffusing in the empty space. Once they finally find an oppositely 



Simulations of the RPM at low temperature and density 



5 



charged ion, they strongly bind to it forming a neutral dimer (or a higher cluster) that 
rarely breaks. As a consequence, the computational time needed to equilibrate the 
system can be astronomically long. To improve the equilibration time in our NVT MC 
scheme, we have adopted two established MC moves and implemented a new one: 

• small and large particle displacements; 

• small and large cluster displacements; 

• formation and breakage of clusters. 

3.1. Small and large particle displacements 

The first move we select is the standard single-particle displacement, where the x, 
y, and z coordinates of a randomly selected particle are each displaced by a small 
amount 5 chosen randomly from the interval {— <5 max , <5max}- The move is accepted with 
the standard Metropolis probability min (1, e~^ u( - n ^~ u ^), where u(o) and u(n) are the 
energies of the particle before and after the trial move, respectively, and j3 = 1/ksT. 
According to normal practice, <5 max can be adjusted to give some desired acceptance rate 
for the move over the course of the simulation. In principle, such moves should allow each 
particle to diffuse as a free ion, and join or leave a cluster. However, when the density 
is very low (and the simulation box is very large), short single-particle displacements 
are not sufficient to sample the phase space properly. Thus, we also randomly attempt 
displacements where <5 max — L/2. These occasional large displacements are intended to 
accelerate cluster formation (if thermodynamically favourable) and to allow the system 
to explore more significant regions of phase space within the simulation timescale. 

3.2. Small and large cluster displacements 

Single-particle MC moves are not enough to equilibrate highly clustered systems, 
and so we also implement a cluster move similar to that proposed by Orkoulas and 
Panagiotoupoulos [11]: we first identify all of the clusters in the system, then choose a 
cluster at random and select displacements from either a small or a large interval, as 
in the single-particle moves. Cluster moves that result in the merging of two or more 
clusters have to be treated extremely carefully in order to respect detailed balance; the 
reverse move has to be attempted with equal probability as the forward move. Here we 
take a simple solution, and simply reject all cluster moves that lead to the merging of 
clusters [11,27]. In this way, the instantaneous cluster distribution is left intact, and 
detailed balance cannot be violated. With this simple approach, the cluster move is 
accepted with the normal Metropolis probability min (1, e~^ u( - n ^~ u ^), where U{6) and 
U(n) are the energies of the system before and after the trial move, respectively. Of 
course, this move does not lead to the formation or breakage of clusters; a specific move 
to effect these transformations is detailed next. 
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3. 3. Novel move for the formation and breakage of clusters 

The last attempt we make to improve the equilibration of the system is to introduce a 
novel move that offers the opportunity of forming and breaking clusters. This 'cluster 
formation/breakage' (CFB) MC move is designed to respect detailed balance, and is 
implemented as follows: 

(i) we choose a particle at random (particle 1), without knowing a priori whether it 
belongs to a cluster; 

(ii) we identify all of its neighbours within a cut-off distance A, which can be tuned to 
give optimal performance, as described below; 

(iii) we choose a neighbour at random (particle 2), irrespective of its charge, and store 
its separation from particle 1, r 12 (o); 

(iv) we then move particle 2 to a new separation from particle 1, r 12 (n), chosen randomly 
and uniformly from the interval a < r 12 (n) < A, and with a random orientation of 
the corresponding separation vector; 

(v) we accept the move with a probability min(l, [ri2(n)/ri2(o)] 2 e~^' 7 ' ri ' )_f/( ' ^), where 
U(o) and U(n) are the energies of the system before and after the trial move, 
respectively 

This CFB move respects detailed balance, and does not add any bias towards the 
formation or breakage of a cluster, since clusters can form and break, and isolated 
ions can simply be displaced (see Appendix). 

In the current work, the target acceptance rate for the single-particle and cluster 
moves is approximately 40%, while the acceptance rate for the CFB moves varies 
between a few percent (at high density) and 40% (at very low density). 

4. Results 

We start by comparing our simulation results (using single-particle, cluster, and CFB 
moves) with those obtained by Allahyarov et al. [36] under the conditions T = 300 K, 
e = 8, q = e, and a = 10 A, corresponding to a reduced temperature T* ~ 0.14. 
The molar concentrations of salt lie in the range 10" 4 -10- 2 mol L" 1 . Allahyarov et al. 
counted "the number of oppositely charged pairs which are closer than 3cx", whereas 
we consider associated ions belonging to clusters with two or more ions (with the same 
cut-off). Figure 1 shows the total concentration of associated ions p a as a function of 
the total ion concentration p. Data from figure 1 of reference [36] are included after 
multiplying the salt concentration and associated ion-pair concentration by two; in 
keeping with these data, we quote the concentrations in units of mol L _1 , assuming a 
particle diameter a = 10 A. All we want to emphasise here is that we get good agreement 
with the established results for concentrations in the range 10 _4 -10~ 2 mol L~ . 

After having confirmed that the algorithm is working in the density regime that has 
already been studied, we move to the central aim of this manuscript, i.e., equilibrating 



Simulations of the RPM at low temperature and density 



7 



0.2 



0.15 ■ 



■- -■ X=3o This work 
o X=3a Ref.[31] 



2 o.i - 



0.05 - 



0.05 



0.1 



0.15 



p (mol L" ) 



0.2 



Figure 1. Concentration of associated ions p a against the total ion concentration p. 
The open circles are results from reference [36] and the red squares are our results. 
The cut-offs chosen to identify the clusters are indicated in the legend. 



the RPM at extremely low density. The lower the temperature and the density, the 
longer it takes to equilibrate the system. Thus, a good test for the algorithm is to 
equilibrate the system at the lowest density and lowest temperature of interest, i.e., 
T* = 0.04 and p* = 9.03 x 10~ n . From preliminary tests, it appears that at the same 
temperature and higher densities of around p* = 10~ 6 the system equilibrates in a 
reasonably short time. 

For clarity, we define three different MC cycles: Monte Carlo cycle (MCO) consists 
of N moves, 90% of which are small displacements of single particles, 3% are large 
displacements of single particles, 3% are are small displacements of randomly chosen 
clusters, and 4% are large displacements of randomly chosen clusters; Monte Carlo 
cycle 1 (MCI) consists of N moves, 90% of which are small displacements of single 
particles, 3% are large displacements of single particles, 3% are small displacements of 
randomly chosen clusters, 2% are large displacements of randomly chosen clusters, and 
2% are CFB moves; Monte Carlo cycle 2 (MC2) consists of N moves, 70% of which are 
small displacements of single particles, 10% are large displacements of single particles, 
5% are small displacements of randomly chosen clusters, 5% are large displacements of 
randomly chosen clusters, and 10% are CFB moves. In all cases, we define a cluster 
according to A = 2a; in MCI we choose A = L/2, whereas in MC2 we consider different 
values of A. 

In figure 2 we show the reduced density of associated ions, p*, versus MC cycle for 
simulations run according to the MCO, MCI, and MC2 protocols, and starting from the 
same initial configuration. MCO shows almost no structural evolution on the simulation 
timescale, and hence is entirely inadequate for simulations at low temperature and 
density. This is caused by the incredibly long distances an ion should cover in order to 
find another ion (difficult cluster formation), and at the same time by the low probability 



Simulations of the RPM at low temperature and density 



8 



of thermally activated dissociation of ion pairs at very low temperature (difficult cluster 
breakage). It is evident that all of the MC2 runs and the MCI run equilibrate to the 
same structure, within our simulated time scale. Moreover, the equilibration times are 
quite different: MC2 (with 10% CFB moves) equilibrates faster than MCI (with 2% 
CFB moves). 
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Figure 2. Reduced density of associated ions versus Monte Carlo cycle for the 
MCO, MCI and MC2 protocols and the same initial configuration at T* = 0.04 and 
p* = 9.03 x 10 -11 . The legend indicates the cut-off A chosen for CFB. 



Next, we select the MC2 Monte Carlo scheme and equilibrate the system using 
different values of A, to show that the final equilibrium state does not depend on the 
choice of A, but that its equilibration rate does. To this end, we use different values of 
A, ranging from 1000a up to half of the box length L/2 (7070a), and plot the reduced 
density of associated ions versus MC cycle. Figure 2 shows that all of the chosen 
values of A lead to the same equilibrium density of associated ions. Strikingly, the 
equilibration rate decreases with increasing A. Choosing a small value for A allows for 
a faster equilibration; however, A cannot be too small compared to the mean separation 
of clusters, as it will lead again to inefficient sampling, not allowing clusters to merge 
or break. Therefore, the optimal value of A should decrease with increasing density. 

We now demonstrate that the convergence of the algorithm does not depend on 
the initial configuration chosen, and that the system is quasi-ergodic on the simulation 
time scale. To this end, we set A = 1000cr and compute the density of associated 
ions in simulations starting from three completely different initial configurations: (a) a 
configuration containing only isolated ions; (b) a configuration containing 40% isolated 
ions and 60% ions in pairs; and (c) a configuration containing 25% isolated ions and 
75% ions in pairs. Figure 3 shows that convergence is achieved irrespective of the initial 
configuration. It is also encouraging that the algorithm allows for significant fluctuations 
in the number of associated ions, which indicates that there is a dynamic equilibrium 
involving the formation and breakage of clusters. 
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Figure 3. Reduced density of associated ions versus Monte Carlo cycle in simulations 
started from configurations (a), (b), and (c) (see text), at T* = 0.04 and p* = 
9.03 x lO -11 . 



5. Conclusions 

In this manuscript we have presented a numerical method for speeding-up computer 
simulations of the restricted primitive model at low temperatures (around T c ) and 
very low reduced densities (down to 10~ 10 ). Our method involves the combination of 
conventional single-particle and cluster moves with a novel 'cluster formation/breakage' 
move, designed specifically to equilibrate the system in a reasonable time, even at 
such extreme thermodynamic conditions. The suggested Monte Carlo scheme is 
straightforward to implement: after having set the value of the maximum neighbour 
distance A the method is inherently efficient, in that the system quickly converges to 
its equilibrium state. This method might also be applied to the equilibration of other 
systems that form aggregates at low concentrations. We should mention that we are 
aware of other techniques that might prove useful to equilibrate very low concentration 
systems, such as the 'geometric cluster algorithm' by Liu and Luijten [43], demonstrated 
to speed up simulations of complex fluids near criticality and/or with differently sized 
components, and a novel cluster move by Almarza [44,45]. As far as we are aware, our 
results extend to far lower concentrations than in any previous studies on the vapour 
phase of the restricted primitive model. The algorithm presented here allows for a 
comprehensive study of the vapour phase around the critical temperature and at reduced 
densities down to 10~ 10 : such low densities seem to be relevant for experiments on low- 
concentration solutions of ions in low-dielectric organic solvents. A detailed report of 
our investigations is in preparation. 
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Appendix: Acceptance move of the cluster formation/breakage move 



Below we derive the acceptance rule for the CFB move, and show that it satisfies the 
detailed balance condition. Detailed balance requires that 

p(o)ii(o — > n) — p{n)n{n — > o) (1) 

where p(o) is the probability that the system is initially in the old configuration o, 
7r(o — > n) is the transition probability from the old to the new configuration n, p(n) 
is the probability the system is initially in the new configuration, and ir(n — > o) is the 
transition probability from the new to the old configuration. Each transition probability 
in equation (1) can be expressed as the product of two terms: 

7r(o — y n) = a(o — > n) x accio — > n). (2) 

a(o — > n) is the probability of generating a new configuration n starting from o, and 
acc(o —> n) is the probability of accepting the move. A similar equation holds for 
ir(n — > o). In our simulations, the old configuration is defined by choosing two particles 
(1 and 2) at random, and computing their relative distance ri 2 (o), and the total energy 
of the system U(o); the new configuration is generated by displacing particle 2 with 
respect to particle 1, and computing their new relative distance ri2(n), and the new 
total energy of the system U(n). r 12 (n) is generated uniformly on the interval {a, A}, 



and hence a(o — > n) 
Combining equations 
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and (2) gives 
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To conclude, we implement a Metropolis sampling scheme using an acceptance 
probability for a move from o to n of 
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